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Abstract 

An Euler solution for an axisymmetric jet engine con- 
figuration without blade effects is presented. The Euler 
equations are solved on a multiblock grid which covers 
a domain including the inlet, bypass duct, core pas- 
sage, nozzle, and the far field surrounding the engine. 
The simulation is verified by considering five theoretical 
properties of the solution. The solution demonstrates 
both multiblock grid generation techniques and a foun- 
dation for a full jet engine throughflow calculation. 

1. Introduction 

Numerical simulations increasingly realize the hope of 
simulating a wide range of physical processes for engi- 
neering design and scientific research. One facet of re- 
alizing this goal is developing numerical simulations for 
geometrically complex domains. Another facet is devel- 
oping numerical techniques appropriate to a problem of 
interest. Thus the simulation of an engine configuration 
is both a testbed for simulations for complex geometries 
and provides a tool for engine research and design. 

A central issue in performing simulations in complex ge- 
ometries is determining the grid: the choice of points on 
which the equations will be discretized and how these 
points are organized to interact with the numerical solu- 
tion techniques. There are many specific approaches to 
this problem, but the grids are either locally structured 
or unstructured. Structured grids have a trivial local re- 
lationship between neighbors while unstructured grids 
have neighbor relationships stored in tables. This dif- 
ference has implications for the algorithms which may 
be implemented, the methods of generating grids and 
ancillary issues. For example, matrix formulations of 
implicit schemes on unstructured grids lose the band- 
edness possible with structured grids. However, when 
using structured grids, the burden of dealing with com- 
plex geometries shifts to how the blocks interact with 
each other. 

For these classes of grids, there are many approaches 
to complex geometry simulations. Euler solutions for 
aircraft configurations using unstructured grids have 
been presented by Jameson, Baker, and Weatherill 1 and 


This paper is declared a work of the U.S, Govern- 
ment and is not subject to copyright protection in the 
United States 


Lohner and Parikh. 2 Buning et. al. 3 have simulated 
an ascending Space Shuttle using overset grids, and 
Chesshire and Henshaw 4 have demonstrated Navier- 
Stokes solutions with overlapping grids in two dimen- 
sions. Non-overlapping, structured multiblock grids and 
solutions for aircraft configurations have been used by 
Sawada and Takanashi 5 while Whitfield et. al. 6 have 
found Euler solutions for counter-rotating propfans whe- 
re the structured grids are in relative motion. Hall et. 
al. 7 have demonstrated solutions for a ducted prop-fan. 
Stewart 8 has used multiblock grids to find Euler solu- 
tions for multi-element airfoil sections. 

However, for multi-stage turbomachinery, significant ge- 
ometrical and physical complications are added to a 
simulation. Blade rows are in relative motion and in- 
teract strongly so simplifications and physical modelling 
are still necessary. One pioneering approach to this 
problem was proposed by Wu 9 , and it represents the 
flow through relatively moving blade rows as two dis- 
tortable streamsurfaces. Jennions and Stowe 10,11 have 
shown quasi-three dimensional throughflow calculations 
for a turbine stage. Wisler, Koch, and Smith 12 mention 
the use of circumferential average flow calculations in 
preliminary compressor design. More recently, averag- 
ing models 13 ' 14,15 have been proposed to simulate three 
dimensional multi-stage turbomachines. 

Although component performance is important, there 
are interactions between compressors, combustors, and 
other components. After initial analysis and experi- 
mental testing of an engine design, components are in- 
tegrated based on their performance in isolation. Thus 
consideration of engine integration effects is limited to 
one dimensional analysis and final engine testing, even 
when less conservative designs require consideration of 
these interaction effects. 

The purpose of the current work is to lay a founda- 
tion for modelling some of these integration effects. As 
a step in this task, the current calculation deals with 
the geometrical and numerical problems of combined 
internal and external flow for an engine configuration. 
Section 2 details the geometry modelling and motivates 
the choices and changes which have been made. Sec- 
tion 3 summarizes the techniques used to create the 
multiblock grids, and Section 4 explains the numerical 
methods employed. The numerical results are presented 
in Section 5. 
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2 . Geomet ry Modelling 

The geometry of the simulated engine configuration is 
based on the General Electric Energy Efficient Engine 
( E 3 ). The geometry is specifically based on the Inte- 
grated Core/Low Spool (ICLS) ground test configura- 
tion for which there is extensive test data. 16 Develop- 
ment of the E 3 engine was sponsored by NASA as part 
of a broad-based effort to increase aircraft efficiency. 

Several modifications have been made to this geometry 
to achieve the desired configuration and to allow for 
a convergent solution. No attempt has been made to 
model the combustor of the engine at this stage of the 
simulation, and a duct has been installed from the com- 
pressor to the turbine section. A small booster island 
has been excluded. Further, because the ICLS ground 
test did not include the external flowpath of the na- 
celle, one has been designed for this case. An ellipse is 
used from the throat of the inlet to the highlight. The 
external nacelle contour is a NACA 1-Series forebody 
and the aft section is an arbitrarily-fit Bezier spline. In 
the future, improved geometrical descriptions may be 
readily incorporated into the simulation. 

The surfaces are modelled using splines parameterized 
by arc length. To ensure the integrity of this represen- 
tation, continuous first derivatives of the raw data are 
ensured wherever possible. Also, it has been necessary 
to smooth rough surface data in the compressor and 
turbine sections. A least squares fit to a basis set of 
sine, cosine functions was successful in removing these 
high frequency geometrical irregularities. 

In the current calculations, blade effects have been ne- 
glected while the geometrical issues of throughflow cal- 
culations are studied. Without blade effects to do work 
and extract energy from the stream, the flow through 
the core section is significantly changed. In particu- 
lar, the compressor blades do no work to assist the flow 
through the core section and achieve the 23:1 pressure 
ratio of the engine’s compressor. It was necessary to 
change the core surface contour so that the maximum 
to minimum area through the core section is approxi- 
mately 3:1. This change was achieved by both lifting 
the entire splitter and altering its core contour. 

3. Grid Generation Techniques 

The grid for this engine configuration is developed using 
the TOPOS program. 8 ’ 17 This program determines a 
domain decomposition such that the domain is covered 
with non-overlapping regions which are each topologi- 
cally rectangular and hence called blocks. Each block 
is dimensioned so that it contains a structured grid 
while the block interfaces have coordinate lines continue 
through them without slope discontinuities. 

The domain decomposition is performed using a search 


algorithm which finds boundary conforming regions in 
a two-dimensional domain. 17 In the same way that the 
skin of a balloon will conform to the bounding walls 
when blown up in a confined space, this algorithm re- 
fines a coarse approximation to the perimeter of a region 
so that it conforms to any nearby neighboring bound- 
aries without excessive stretching. 

This algorithm is demonstrated in figure 1 by consid- 
ering the simple case of an airfoil in a box, PQRS , 
and an initial, coarse approximation to the region be- 
low the airfoil, ACRS. To transform ACRS to a region 
which conforms to the lower surface of the airfoil, the 
algorithm finds the point which determines the high- 
est flat ceiling above the middle third of AB , EF, as 
in figure lb). Limiting the depth of this probe to the 
height Search Depth prevents finding points on the seg- 
ment PQ which would yield a perimeter with excessive 
stretching. The probe considers the curves which de- 
fine the domain and finds the point D. The perimeter 
is then modified to ADBC, and two child segments, 
AD and DB, are created from AB. 
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Figure 1: Steps in Finding a Boundary Conforming Re- 
gion. 

The refinement of ADB involves refining AD and DB. 
AD is refined, but DB is not since it is part of a bound- 
ary defining curve. Probing above the middle third of 
AD, GH , yields no point within search depth so the 
two outer segments AG and HD are considered for re- 
finement. This failure to find a highest ceiling implies 
there will be a transition between bodies in the perime- 
ter. This refinement is extended in a tree structure by 
considering successively smaller sub-segments as above, 
and eventually yields the perimeter AZB of figure Id). 
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With this algorithm, boundary conforming regions are 
successively found and removed until the domain is cov- 
ered. 

After finding a decomposition, transformations may be 
applied to it which involve cutting and merging blocks 
in general ways. The data structures which describe 
the decomposition are designed in such a way that these 
transformations are isomorphic with respect to the data 
structures, providing a tool for the development and 
adaptation of multiblock grids. 

Given a domain decomposition, it is necessary to de- 
termine the dimensions of the grids in each block. To 
have a structured grid with coordinate lines continuous 
through the block interfaces, two conditions must be 
satisfied. First, the numbers of cells on opposite sides 
of a block must be equal; second, the number of cells 
on opposite sides of a block interface must be equal. 
These conditions may be expressed as an undercon- 
strained system of linear, integer coefficient equations 
with integer solutions. If the equations have an admis- 
sible solution, a linear programming algorithm will find 
the solution. Further, the free parameters in the so- 
lution are also found, providing a technique for global 
adaptation of the grid. 

To obtain grids within the blocks, elliptic grid genera- 
tion 18 is used to smooth initial approximations to the 
grid found by algebraic methods. The interfaces be- 
tween blocks are also smoothed. Having a locally smo- 
oth grid which may be thought of as a locally C 1 trans- 
formation between physical and computational space 
results in a second order local discretization error in 
the numerical approximation, as will be explained in 
Section 4. 


4. Numerical Methods 


The two-dimensional Euler equations model inviscid, 
compressible, rotational flow and allow entropy and vor- 
ticity production across shocks. For future applications, 
they will allow for body forces and energy addition. The 
Euler equations are given by 


dt 


J wdV ol = 
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Here 7 = 1.4. The circumferential velocity is assumed 
to be identically zero. 


The numerical approximation and solution techniques 
are based on the FL052 program of Jameson et. al . 19 
The equations are discretized by taking the dependent 
variables, w , to be located at the centroids of the grid 
cells. The convective flux for a cell, Q(w ), is approx- 
imated from the average dependent variable value be- 
tween adjacent cells and the area of the cell face. 

The discretized equations for each cell are advanced to 
a steady state by a multi-stage scheme: 


= Wi + ai A t(Q(wi) + D(wi)) 
w ^ =Wi + £*2 At(Q(w^) + D(w^)) 
w = Wi + CK3A t(Q(w^) + (2) 

^(4) = w . _j_ a4 /\t(Q(w + D(w^)) 
wi + 1 = Wi + a 5 At(Q(w {4) ) + D(w (1) )) 


where a* are coefficients, Q(w) is the convective flux ap- 
proximation for the cell, and D(w) is the artificial dissi- 
pation. The artificial dissipation consists of third-order 
dissipation, which stabilizes the time stepping scheme, 
and first-order dissipation, which is switched on near 
shocks to capture them. 

Neumann stability analysis of Equation ( 2 ) for a one- 
dimensional analog of Equation ( 1 ) indicates that sta- 
bility is expected for appropriate coefficients cxi. Since 
a steady state solution is sought, convergence is acceler- 
ated by using the maximum time step in each cell. For 
isenthalpic conditions, convergence may be enhanced 
by using enthalpy damping . 20 Enthalpy damping uses 
the Total Energy, E> and kinetic energy, ~(u 2 + v 2 ) } 
to estimate the local deviation from constant enthalpy 
and force the mass and momentum equations accord- 
ingly. In the absence of body forces and energy release, 
enthalpy damping has been used. 

Convergence has also been enhanced by using a multi- 
grid algorithm . 21 Inset grids are included by having 
each grid dimension contain a factor of a power of two. 
The number of multigrid levels is limited by the coars- 
est unstretched grid which in turn is determined by the 
shortest block side. For this calculation, three levels can 
be constructed, and their use improves convergence con- 
siderably. The rate of convergence, measured by the av- 
erage log residual reduction per cycle, doubles for each 
added level. 


f{w) = 
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At solid boundaries, a free slip condition is applied and 
pressure is extrapolated to the surface by considering 
grid skew and centrifugal effects. At the boundary in 
the far field, a uniform free stream velocity is assumed 
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Figure 2: Inner Region of the Grid for the Jet Engine Configuration 


and matched to the solution using Riemann invariants. 
At interfaces between blocks, a simple matching of the 
solutions is possible. Since coordinate lines are contin- 
uous through the interface, the grid may be continued 
across a block interface. Consequently, the dependent 
variable values are known for the flux approximation, 
Q(w), and artificial dissipation, D(w) at each interface. 
Because communication occurs between grids by copy- 
ing data values at each stage of Equation (2), interpo- 
lation is unnecessary. Further, these conditions result 
in the property that splitting blocks does not change 
the solution or convergence — they are identical in both 
cases. The only penalties paid for splitting blocks are 
for setting up additional loops and subroutine calls and 
for decreased vector lengths. 

Since coordinate lines are C 1 through an interface, it 
is possible to make theoretical arguments that the so- 
lution is locally second order accurate. 17 The trunca- 
tion error of the flux approximation, <3(xe), may be ex- 
pressed on the locally uniform computational domain 
because there is a local C 1 transformation between phys- 
ical and computational space. In this form, it can be 
shown that the truncation error should be second or- 
der. This result has been numerically verified on grids 
of this type. This conclusion is also supported by ex- 
amination of the drag coefficient computed on grids of 
varying resolution. 


5, Results 

The geometrical domain includes the contours for the 
hub, splitter, and nacelle which are symmetric about 
the engine axis, and extends 13 engine lengths into the 
far field to a circular boundary of truncation. The grid 
for the domain, shown in figure 2, contains 16800 cells 
in 35 blocks whose edges are indicated by bold lines. 

The grid contains three inset multigrid levels which are 
used to enhance convergence. These grid levels are cre- 
ated from a coarse base grid by multiplying each grid 
dimension by a factor of four. Additional grid levels 
would enhance convergence further, but using a coarser 
base grid with this non-uniform grid topology would 
result in excessive grid stretching and degrade the ac- 
curacy of the solution. 

Although the geometry is symmetric about the center- 
line, the grid is not. The grids in the upper and lower 
half planes are generated without a symmetry condi- 
tion applied between the two halves, except that the 
engine axis corresponds to a coordinate line. This grid 
asymmetry is exploited to verify the solution, since the 
solution will not be symmetric because of the grid but 
because of proper resolution. Because the coordinate 
line on the axis does not necessarily correspond to a 
block interface, some blocks are in both the upper and 
lower half planes. Having the centerline correspond to 
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Number of Cycles 


Figure 3: Convergence Histories with Two Multigrid 
Levels (solid line) and Without Multigrid (dashed line) 

the engine axis results in a condition of zero flux through 
the axis, which isolates the two half planes. There is 
also circumferential isolation since the grid does not ex- 
tend in the third dimension, theta. 

The Euler solution was run for 2500 cycles on this grid 
at Mach M = 0.3 and an angle of attack a = 0°. The 
Mach number, M = 0.3, was chosen for this case to 
avoid transonic flow and choking in the engine core 
passage. Since no work is done by the blades on the 
fluid through the compressor section, the fluid accel- 
erates dramatically and can become sonic only at the 
throat of the duct. This sonic point limits the mass flow 
through the duct and hence the capture streamline for 
the core. A low Mach number was chosen to minimize 
this problem. 

Although the ground test data 16 is not an appropri- 
ate basis for comparison at this stage because no blade 
effects are being considered, the results may be veri- 
fied based on five theoretical properties of the solution. 
First, the solution is convergent in the sense that mass 
is conserved locally. The density residual, RMS av- 
eraged over the field, is shown in figure 3, and indicates 
convergence. This quantity is a direct measure of how 
close to mass conservation the solution is. The energy 
equation is also convergent. Convergence of the total 
energy residual, over the field is a direct measure 
of how close to energy conservation the solution is. 

Second, since the solution is subcritical, it should be 
isentropic, and deviations are an indication of error. 
The deviation from constant entropy over the surfaces 
is shown in figure 4. The leading edge of the nacelle 
contributes the largest amount to the entropy devia- 
tions in part because the limited ingestion of air into 
the engine inlet moves the capture streamline closer to 
the engine axis and the stagnation point into the inlet. 
As a result, substantial turning of the flow around the 
leading edge of the nacelle is required. The leading edge 



Figure 4: Fractional Deviation of Entropy From its Free 
Stream Value Over the Surfaces of the Configuration 


5 





Pressure Coefficient Pressure Coefficient Preenue Coefficient 

1.20 0.80 0.40 0.00 - 0.40 - 0.80 - 1.20 - 1.60 1.20 0.80 0.40 0.00 - 0.40 - 0.80 - 1.20 - 1.60 1.20 0.80 0.40 0.00 - 0.40 - 0.80 - 1.20 - 1.60 



* Inner Surf see 
. Outer Surface , 

l 



Figure 5: Pressure Distributions Over the Symmetric 
Surfaces of the Configuration 


of the splitter plate is difficult to resolve because it has 
a small leading-edge radius, which creates an entropy 
deviation. 

Third, by geometrical symmetry, the lift coefficient, Cu 
should theoretically be zero, and on a symmetric grid 
in symmetric conditions, one would expect the numer- 
ical result to be identically zero due to cancellation. 
However, for a symmetric geometry and an asymmetric 
grid, one cannot expect an identically zero value. The 
computed lift coefficient is C\ = -0.0009 and is due to 
grid asymmetry, the truncation errors in the numerical 
approximations, and artificial dissipation. 

Fourth, since the Euler solution is sub-critical, the drag 
coefficient, Cd, theoretically should be zero. The value 
Cd = 0.0009 can again be attributed to grid asymmetry, 
errors due to numerical approximations, and artificial 
dissipation. 

The symmetry of the geometry and the asymmetry of 
the grid provide a fifth test. The surface distributions 
of flow quantities should be symmetric. The surface 
pressure distributions are sensitive to errors and the 
distributions for the three symmetric components are 
shown in figure 5. 

7. Conclusions 

A solution to the Euler equations for an axisymmetric 
jet engine configuration without blade effects has been 
presented. The solution has been found on a multi- 
block grid which includes both the internal ducts and 
the external flow surrounding the engine. Further, the 
result has been verified on the basis of five theoretical 
properties of the solution. Thus this work demonstrates 
techniques for finding numerical solutions in complex 
geometries and provides a foundation for more detailed 
modelling of engine components. 
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